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The influence of the choice of pharmacokinetic model on subsequent determination of covariate 
relationships in population pharmacokinetic analysis was studied using both simulated and real 
data sets. Simulations and data analysis were both performed with the program NONMEM. Data 
were simulated using a two-compartment model, but at late sample times, so that preferential 
selection of the two-compartment model should have been impossible. A simple categorical covari- 
ate acting on clearance was included. Initially, on the basis of a difference in the objective function 
values, the two-compartment model was selected over the one-compartment model. Only when the 
complexity of the one-compartment model was increased in terms of the covariate and statistical 
models was the difference in objective function values of the two structural models negligible. For 
two real data sets, with which the two-compartment model was not selected preferentially, more 
complex covariate relationships were supported with the one-compartment model than with the 
two-compartment model. Thus, the choice of structural model can be affected as much by the 
covariate model as can the choice of covariate model be affected by the structural model; the two 
choices are interestingly intertwined. A suggestion on how to proceed when building population 
pharmacokinetic models is given. 



KEY WORDS: population pharmacokinetics; structural model ; covariate model; NONMEM ; 
model selection. 



INTRODUCTION 

The population approach for pharmacokinetic data analysis, as imple- 
mented within the program NONMEM (1), was originally proposed for the 
analysis of sparse data (2). As such, it was not originally intended to be 
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used as a tool for selection of the structural (pharmacokinetic) model. The 
choice of structural model was to be based upon the best a priori information 
about the drug to be studied, or a simple model was to be used when the 
nature of the data was such that a more complex model could not be sup- 
ported. Population pharmacokinetic analysis, having gained widespread 
acceptance, is used for analysis of both sparse and more intensively sampled 
data sets, but the initial selection of the structural model is still somewhat 
arbitrary. A review of 30 pharmacokinetic data analyses (3-32) in which 
NONMEM was used revealed that only one-third reported statistical just- 
ification for the choice of structural model. The investigators that did report 
a reason for their choice of structural model did not consider that the choice 
made may have had any influence on subsequent statistical or covariate 
model building. Whether or not an attempt is made to justify the selection 
of the structural model, there usually exists an implied and untested assump- 
tion that the relationships between "model-independent" parameters (such 
as clearance) and covariates are unaffected by the structural model. 

Simple pharmacokinetic models may be used for a variety of reasons, 
the most obvious being that a more complicated model cannot be justified 
statistically. For example, the limit of assay sensitivity may preclude charac- 
terization of longer terminal phases of drug elimination, or data obtained 
at steady state may be insufficient to enable characterization of a distribution 
phase. Within the field of therapeutic drug monitoring, simple pharmaco- 
kinetic models are routinely used for prediction purposes. For example, 
samples taken for the optimization of digoxin therapy are intentionally 
obtained after the distribution phase. Again, in these different situations, 
the covariate-pharmacokinetic relationships are assumed to be unaffected 
by the pharmacokinetic model itself. 

This paper seeks to examine the influence of the choice of structural 
model on the subsequent building of the covariate model. Initially, we 
studied simulated data sets (see "Simulations"). Analysis of the simulated 
data revealed that the issue was more complex than originally thought and 
led us to reconsider criteria for choosing the structural model itself. Analysis 
of two real data sets (see "Real Data") exemplify the issues encountered 
with the simulations. Finally, the implications of the results are discussed. 

SIMULATIONS 
Data 

Data were simulated using the program NONMEM (Version IV), using 
a two-compartment model with intravenous input (ADVAN3, TRANS3) 
(I). The values of the pharmacokinetic parameters were loosely based upon 
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those of digoxin. For each of 100 subjects, clearance (CL) was sampled from 
a normal distribution characterized by a population mean of 11.3 L/hr, 
and constant coefficient of variation ((o C i) of 25% ; volume of the central 
compartment (Vc) was sampled from a normal distribution characterized 
by a population mean of 37.6 L, and (o Vc of 30%. No covariance between 
CL and Vc was included. Intercompartmental clearance (Q) and volume of 
distribution at steady state (Vss) were fixed to be 31.4 L/hr and 470 L, 
respectively, for all subjects. The values of the pharmacokinetic parameters 
were such that the typical values for /| /2 < a ) and U nm were 35 min and 36 hr, 
respectively. A proportional random normal error (10%) was added to each 
concentration. 

A second series of data was also simulated with a simple dichotomous 
covariate (e.g., presence of an interacting drug) affecting CL in a subset of 
individuals within each data set. Two subpopulations were created, one (n = 
75) with a mean CL of 1 1.3 L/hr (e.g., monotherapy), and another (w = 25) 
with a mean CL of 14.2 L/hr (e.g., in the presence of an interacting drug). 

The goal was to simulate data after a borderline time point, with which 
a two-compartment model could no longer be preferentially selected. Thus, 
data within four different time ranges (3-120 hr, 6-120 hr, 12-120 hr and 
48-120 hr) were tried. Each subject had three samples selected randomly 
from a series of fixed time points within the time ranges. Fifty data sets, 
each consisting of 100 subjects, were generated for each time range and for 
each of the two situations, without and with the drug interaction covariate 
influencing CL (in a subset of individuals). 

Analysis 

Data analysis was carried out in two stages. First, to find the borderline 
time point (after which a two-compartment model could no longer be just- 
ified statistically with the simulated data), one- and two-compartment mod- 
els were fit to all the data sets. The analysis included a covariate effect only 
when the effect was simulated in the data. Second, after the borderline time 
point was identified, different covariate models with both the one- and two- 
compartment models were fit to the data simulated after this point. 

Model selection was based on the objective function value of the reduced 
model minus that of the full model, which approximates the difference in 
minus twice the maximum log likelihood of the two models. The difference 
is approximately chi-square distributed. 

Enigma 

The frequency of selecting the two-compartment model, with and 
without the covariate in CL, was much larger than would be expected for 
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Lower bound of time range simulated (h) 

Fig. 1. The percentage of data sets simulated between the time indicated and 120 hr, in which 
a two-compartment model was selected (based on p^0.05) over a one-compartment model. 
The squares are the percentages associated with the data that included the drug interaction 
covariate on CL in both the simulation and analysis; the triangles are the percentages when 
no covariate was present. 

each of the four time ranges studied (Fig. 1). For example, because the 
t\/2(a) was only 35 min, data simulated between 48-120 hr were over 80 half- 
lives after the end of the distribution phase, making the selection of the two- 
compartment model appear nonsensical. No matter what the borderline 
time point was, the two-compartment model was virtually always selected, 
making the objective of finding the borderline between the one- and two- 
compartment models impossible to attain. 

Resolution of the Enigma 

For each set of two-compartment parameter values used in the simula- 
tions, we determined the values of CL and volume ( V) corresponding to 
the exponential associated with the terminal slope (analogous to treating 
data originating from the terminal phase of a two-compartment model as if 
it originated from a one-compartment model). Two findings resulted from 
these calculations. First* the drug interaction covariate influenced V as well 
as CL. Second, a high correlation between CL and V was found. (Recall, 
the simulation of two-compartment parameters assumed no correlations.) 
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Table I. Mean Difference (and SD) in the Objective Function Values of the One- and Two- 
Compartment Models, for the Fifty Replicate Data Sets" 



True covariate Covariatc Statistical 

model model model 

(Simulation) (Analysis) (Analysis) 



Time range of simulated data (hr) 



3-120 6 120 



12-120 48-120 



No covariate No covariate No covariance between CL 

onCX on CL and ft- (f) 272(41) 99.5(18) 60.1 (19) 7.1 (25) 

No covariate No covariate Covariance between CL 

on CL onCL and Vc (V) 235(39) 10.9(7.7) 0.6(2) 0 (0) 

Covariate on Covariate on No covariance between CL 

CL h CL andKr(K) 300(41) 112 (24) 71.7(23) 19.3(28) 

Covariate on Covariate on Covariance between CL and 

CL CL Vc (V) 263(40) 24.4(16) 16.4(12) 5.0(5.8) 

Covariate on Covariate on Covariance between CL 

CL CL and V r and Vc (V) 250(40) 8.6(11) 1.3(5.7) 0.6(3) 

•Regarding the difference in degrees of freedom between the one- and two-compartment models as two, a 
difference in the objective function values of 5.99 corresponds to a p of 0.05. 

*For data sets simulated with a simple categorical covariate in CL* 25% of the subjects had a 25% increase 
m CL. 

Tor the one-compartment model, the covariate acts on CL and V; for the two-compartment model, the 
covariate acts on CL only. 

The data simulated for the purpose of finding the borderline time point 
were then reanalyzed. First, a covariance between CL and V was included 
in the one-compartment model. In the case where no covariate had been 
included in the simulated data, inclusion of the covariance between CL and 
V was sufficient to enable the one-compartment model to be selected (Table 
I). With the simulated data that contained the drug interaction covariate, 
the two-compartment model was still often selected, although now at a 
significantly reduced frequency. Second, in the case where the drug inter- 
action covariate was simulated, the covariate was also introduced into V of 
the one-compartment model. Inclusion of the covariate on V finally brought 
the frequency of choosing the one-compartment model closer to that reason- 
ably expected. 
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Results 

Based on the results presented in Table I, data simulated with the 
covariate on CL and within the time range 12-120 hr were chosen to investi- 
gate covariate detection. Regardless of the covariance structure, the fre- 
quency of detecting the covariate effect, when it truly existed, was always 
higher using a two-compartment model than a one-compartment model 
(Tables II and III). When no covariance between CL and Vc (V) was 
included, the rate of detecting the covariate on CL was 87% for the two- 
compartment model and 72% for the one-compartment model. When a 
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Table II. Percentage of Two-Compartment Runs in which the Criterion for Detecting the 

Covariate was Satisfied" 



% of runs selecting a covariate* 


Statistical model 


on 


on on Vc in presence 


on 


on Km in presence 


(Analysis) 


CL only 


Vc only of covariate on CL 


Vss only 


of covariate on CL 


No covariance between 










CL and Vc 


87 


15 4 


2 


4 


Covariance between 










CL and Vc 


83 


52 8 


2 


5 



'All results are based on a I degree of freedom chi-square test; a significance level of 0.001, corresponding 
to a difference in objective function values of 10.83, was used. 

*50 data sets were simulated, in each of which 25 of the 100 subjects had a 25% higher CL. 



covariance between CL and Vc (V) was included, the rate of detecting the 
covariate on CL was 83% for the two-compartment model and 58% for the 
one-compartment model. 

The mean (±SE) estimates of CL in the two subpopulations (without and 
with the covariate effect) over the 50 data sets, using the two-compartment 
model, were 11. 3 ±0.6 and 14.5 ± 1.4 L/hr, respectively. (Recall, the true 
values are 1 1.3 and 14.2 L/hr, respectively.) These values with the one-com- 
partment model were, however, overestimated, 1 5.8 ± 1 .0 and 1 9.5 ± 1 .4 L/hr. 

Tables II and III also show the rate of falsely detecting the drug inter- 
action covariate as acting on the volume parameters ( Vc and Vss in the two- 
compartment model, and V in the one-compartment model). For the two- 
compartment model, there was a high incidence of false positive covariate 
effects on Vc, but only when the covariate effect on CL was omitted from 
the analysis. When the covariate effect was included on CL, there was a 
considerable difference between the estimates of Vc: 75.2 ±87.0 L (mono- 
therapy group), and 44.0 ± 56.0 L (drug interaction group). Consistent with 
the low rate of false positives, this difference was, however, imprecisely 
estimated as indicated by the large standard errors. The estimates of Vss for 



Table III. Percentage of One-Compartment Runs in which the Criterion for 
Detecting the Covariate was Satisfied 0 



% of data sets detecting a covariate* 


Statistical model 






on V in presence of covariate 


(Analysis) 


on CL only 


on Konly 


on CL 


No covariance between 








CL and V 


72 


10 


66 


Covariance between 








CL and V 


58 


2 


62 



'All results are based on a one degree of freedom chi-square test ; a significance level of 
0.001, corresponding to a difference in objective function values of 10.83, was used. 
6 50 data sets were simulated, in each of which 25 of the 100 subjects had a 25% higher CL. 
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the two subpopulations were very similar, regardless of whether the covariate 
effect on CL was included in the analysis, and is reflected in the low rate of 
false positives. 

With the one-compartment model, there was a high incidence of falsely 
detecting the covariate acting on V, especially when the covariate effect on 
CL was included in the analysis (Table III). Consistent with this finding, 
the estimates of V for the two subpopulations differed, and were precisely 
estimated: 816±16 L (monotherapy group), and 923±34 L (drug inter- 
action group). 



REAL DATA 

To investigate whether the issues raised with the simulated data might 
be reflected in real data, two real data sets, which have been reported previ- 
ously (25,6), were also analyzed. 

Quintdine 

Briefly, 391 serum quinidine concentrations from 139 adult hospitalized 
men, collected for routine clinical purposes, were available for analysis (25). 
The following demographic data were recorded: age, height, weight, race, 
smoking background, alcohol history, a-l-acid glycoprotein concentration 
(a x AGP\ presence/absence of congestive heart failure (CHF), presence/ 
absence of dialysis or hemoperfusion, and high/low creatinine clearance 
(CrCL). Verme et al found that the fit of the two-compartment model was 
no better than that of the one-compartment model (personal communica- 
tion). We confirmed this finding by comparing the two models in the absence 
of covariate relationships, as well as in the presence of the covariates found 
by us to be statistically significant with the one-compartment model. These 
final covariates on CL included minor (difference in objective function values 
of the reduced and full models of 4-10) influences of race, height,.weight, 
CrCL status and dialysis status, and major (objective function difference 
>10) influences of current drinking status and a\AGP. A major influence 
of a \A GP was also found on K. Consistent with the original report, we 
found that a hyperbolic relationship between CL and a t AGP, and between 
Kand a x AGP, could be supported. 

The covariates found to be significant with the one-compartment model 
were then used with a two-compartment model and tested by stepwise dele- 
tion. Height and CrCL status, which were minor influences on CL in the 
one-compartment model, were negligible and major, respectively, using a 
two-compartment model. Most interesting, a x AGP, identified as a major 
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influence on Kin the one-compartment model, was found to have a negligible 
effect on Vc in the two-compartment model. 

All models included covariances among ka, Vc ( V), and CL. There was 
nothing remarkable about the estimates of these covariances, except that 
with both one- and two-compartment models, ka and Vc (V) were highly 
correlated. 

Netilmicin 

The netilmicin data used in our analysis was a combination of the 
primary and validation data sets collected from 102 neonates of 27 to 42 
weeks gestational age (6). Netilmicin was given intravenously or intramuscu- 
larly in a dosage of 3 to 5 mg/kg once daily. Two to 12 netilmicin serum 
concentrations per neonate (344 total) were measured. Most neonates had 
measurements collected at target times of 2, 8, and 16 hr postdose. The 
following demographic data were recorded : gestational age, weight, length, 
age at the start of therapy, age on the day of concentration measurement, 
and creatinine clearance. 

A two-compartment model could be justified statistically to describe 
the data (objective function value decreased 151, relative to that for the one-^ 
compartment model). As with the simulations, the aim of the present analysis 
was to compare covariate models when the choice of structural model is 
equivocal. Therefore, the data set was modified by removing all concentra- 
tions taken before 5 hr postdose, resulting in 217 concentrations from 101 
neonates. The one-compartment model no longer differed statistically from 
the two-compartment model. 

The major covariates affecting CL (objective function difference >30) 
using the one-compartment model, included weight, age on the day of con- 
centration measurement, and gestational age. For K, a major influence of 
weight, and a minor (objective function difference 4-10) influence of gesta- 
tional age were found. As described previously for quinidine, the same cova- 
riate model was then applied to a two-compartment model, and stepwise 
deletion of covariates was carried out. When a two-compartment model was 
used, the influences of both weight and gestational age on Vc were no longer 
found to be statistically significant. 

All models used covariances between Vc (V) and CL. With both the 
one- and two-compartment models, CL and Vc(V) were highly correlated. 

DISCUSSION 

Selection of the "correct" model has usually been based upon (i) statist- 
ical criteria that assess the improvement in fit as the model increases in 
complexity, and (ii) a principle of parsimony. A simpler model is chosen 
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over a more complex model only if statistically justifiable, i.e., if the fit 
obtained with each model is essentially the same. In this context, what is 
usually meant by simpler and more complex models is that the parameter 
space of the simpler model is embeddable in the parameter space of the 
more complex model. Often, one model for covariates can be regarded as 
simpler than another, given the same structural model. A one-compartment 
model is generally a simpler model than a two-compartment model. 

Which overall model, including structural, covariate, and statistical sub- 
models, is correct? Our results show that the correct covariate and statistical 
submodels can differ for two different structural submodels. Moreover, the 
least parsimonious covariate and statistical submodels can be associated 
with the most parsimonious structural submodel. Our calculations and simu- 
lations clearly illustrate and explain how choice of the simpler structural 
submodel can be justifiable statistically, but incorrect, insofar as with this 
choice, both the covariate and statistical submodels contain artifactual 
elements. 

Our analyses of the quinidine and netilmicin data suggest that these 
considerations might also apply to real data. But they also suggest that the 
data analyst is faced with a fundamental uncertainty regarding which model 
is correct. With both analyses, the number of important covariates was fewer 
with the two-compartment model. In particular, with quinidine and the two- 
compartment model, the influence of a x AGP on Vc was negligible. More- 
over, quinidine is known to follow two-compartment pharmacokinetics. 
Therefore, it is possible that when the data are analyzed with a one-compart- 
ment model, the apparent significant influence of <i\AGP on f is artifactual, 
and the use of this model is incorrect. 

Alternatively, with a highly protein-bound drug like quinidine, one 
might expect an influence of a\AGP on volume (33). Under the two- 
compartment model, one might expect an influence on Vss, if not on Vc. 
However, as with Vc, the influence of a t AGP on Vss appears not to be 
small, but it is not significant, which is to say that a considerable influence 
might exist, but the data do not support it. The uncertainty in the magnitude 
of the influence is considerable, and may arise simply because the two- 
compartment model is overparameterized. This situation is somewhat similar 
to the one with our simulations, where the influence of the interacting drug 
on either of the volume parameters for the two-compartment model truly did 
not exist. When we included the influence of the interacting drug covariate on 
the volume parameters of the two-compartment model, the influence on Vc 
was considerable, although not significant, while the influence on Vss was 
very small. In summary, it is also possible that the use of a one- 
compartment model is adequate and that for the most part, the influence of 
a x AGP on V \% "real." 
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Interestingly, if there is a real effect of <Z\AGP on Vss, one would expect 
it to be characterized by a hyperbolic relationship, and with a drug like 
quinidine which has a large volume of distribution, also expect there to be 
no intercept term. This relationship is just the type we observed between 
a { AGP and Vss (two-compartment model) and also between a x AGP and V 
(one-compartment model), lending some credence to the correctness of the 
simpler, one-compartment model. 

Nevertheless, we have demonstrated that, theoretically, there must be 
some degree of artifactual influence of a t AGP on V. So it is also possible 
that the influence of a } AGP on Vis partly artifactual and partly real. If this 
combination of influences were true, it could help explain why the effect of 
a x AGP appeared to be stronger on V (one-compartment model) than it was 
on Vss (two-compartment model). We have not found a report in the litera- 
ture that describes the magnitude of the effect of a x AGP on quinidine's Vss, 
The situation is similar with netilmicin, another drug known to follow two- 
compartment pharmacokinetics. The effect of weight on V may be artifac- 
tual, real, or a combination. There does not seem to be a way to resolve the 
uncertainty. 

Consider the case in which a covariate was found to significantly influ- 
ence more than one parameter (e.g., CL and V) of the one-compartment 
model, but, based on a priori information, its influence on one of these 
parameters (e.g., V) is entirely unreasonable. Then, when a two-compart- 
ment model was fit to these data, no unreasonable significant influence of a 
covariate on a parameter (e.g., Vc) was found to exist. One might com- 
fortably conclude that the one-compartment model was incorrect, and its 
use led to an incorrectly identified covariate. Unfortunately, the ability to 
make such an unequivocal decision was not possible with either real data 
set we analyzed. We are only in a position to caution readers about the 
matter of model choice. Either a one-compartment model with associated 
covariates and covariances will be unquestionably incorrect, or there will be 
uncertainty about their correctness. 

Based on our experience, the following scheme could be useful when 
building models during population analysis. First, a comparision of struc- 
tural models without covariates should be made. The simplest supportable 
structural model should be used for building the covariate model, simply 
because the NONMEM program will run much faster during covariate selec- 
tion. Second, possible influences of covariates on all pharmacokinetic param- 
eters of the simpler structural model should be assessed, because, as we have 
emphasized, covariate effects (whether real or artifactual) can influence more 
than one parameter. Third, if a more complicated structural model seems 
to be a priori equally likely, then the covariate model built with the simpler 
structural model (including marginally significant covariates) should be used 
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with the more complicated structural model. Covariates rejected with the 
simpler model may then need to be reconsidered. Finally, covariates can be 
retested by stepwise deletion. 

This paper has addressed model choice that was made ambiguous by 
censoring the data in the distribution phase. Whether or not the results can 
be extrapolated to models and data types other than the few studied here, 
they do indicate that careful consideration of the structural model must be 
made when building a population model. Optimally, the structural model 
should be defined prior to population analysis, based on experimental data 
(e.g., in Phase I or early Phase II). Further, as our investigation has demon- 
strated, ignoring likely candidates for the structural model, in the interest of 
simplicity, can be consequential in the selection of covariate and covariance 
models. Finally, if prior information about the structural model is weak or 
lacking, then inspection of the interaction between structural model, covari- 
ates, and covariances may hold important information about the appropri- 
ateness of different structural models. 
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